function [ GDP_deflator_m,RGDP_hat_m,RGDP_hat_f ] = GDP_deflatorCES_fun_approx(P_hat_mnj,...
    pi_M_f,pi_l_f,wage_hat,pi_M,pi_l,X_hat_mnj,...
    VA_n_0,VA_mnj_0,VA_hat_n,X_hat_mnjf,VA_fn_0,VA_hat_f)

%% Preamble 

% Author: Chris Evans (UPF)

% Purpose: Creating the real GDP using the GDP deflator instead of CPI
% inflation at the firm and sector level. Then using VA to isolate the GDP
% deflator. 
%
% Following master notes2 page 11,12,13 : This codes takes parts of
% PriceLoop_fun.m and Tradeshare_fun.m to calculate the Real GDP using the
% real GDP deflator instead of CPI.

% Date: 12/10/2018


global alphaT rhoT lambdaT etaT shockT rho sigma sigmaT lambda eta ...
    tol maxit D_hat_n Xi_hat_mnjf Xi_hat_mnj a_hat_f a_hat ...
    varphi sL_n sPi_n sD_n N J FI_J france oneminus_alpha_vector ...
    labour_share_PWT countrysecD firmsecD_sorted start start_sorted ...
    sum_by_country_dummy ISO alpha_WIOD_j alpha_WIOD_francej sector_algorithm 



%% French firm calculation 

% Initialising a matrix we will use to sum over k countries for France
% This is done for the intermediate price index P^M
log_P_hat_mnj=log(P_hat_mnj(:,france)');
log_P_hat_mnjf0_repmat = kron(log_P_hat_mnj,ones(FI_J,1));
log_P_hat_f_weta = pi_M_f.*log_P_hat_mnjf0_repmat;
log_P_M_f_hat = sum(log_P_hat_f_weta,2);
P_M_f_hat = exp(log_P_M_f_hat);
   
log_P_M_f_hat_wlambda = log_P_M_f_hat.*(1-pi_l_f);
% Taking wage to power of 1-lambda and multiplying by labor share
log_w_hat_f = log(wage_hat(france));
log_w_hat_f_wlambda = pi_l_f.*repmat(log_w_hat_f,[FI_J 1]);

% Create total cost: F*1 matrix
log_b_hat_f = (log_w_hat_f_wlambda + log_P_M_f_hat_wlambda);
b_hat_f = exp(log_b_hat_f);



% Multiplying the the expression by a technology (a) shock 

p_hat_mnjf_france = b_hat_f.*a_hat_f;
            
% p_hat_mnj_france_rempat: Replicating the matrix from mjf to mnjf, hence
% expanding it by N countries
p_hat_mnjf_france_rempat = repmat(p_hat_mnjf_france, [1 N]);

% Calculation of Q hat, exports from country m to country n sector j
Q_hat_mnjf_france = X_hat_mnjf./p_hat_mnjf_france_rempat;

% omega_tilde_mnj_france: Weight using initial period value added
omega_tilde_mnjf_france = VA_fn_0./VA_n_0(france);

% RGDP_hat_mnj_france:  Real GDP change for french firms (sum of all french
% firm weighted real GDP hat is the real GDP hat for France) 
RGDP_hat_mnjf_france = omega_tilde_mnjf_france.*Q_hat_mnjf_france;
 
% RGDP_hat_m_france: Summing over all french firm real GDP hat to get the
% real gdp hat for France
RGDP_hat_m_france = nansum(RGDP_hat_mnjf_france(:));

% Clearing to save memory 
clear log_P_hat_mnj  log_P_hat_mnjf0_repmat  log_P_hat_f_weta log_P_M_f_hat P_M_f_hat log_P_M_f_hat_wlambda;
clear log_w_hat_f log_w_hat_f_wlambda log_b_hat_f b_hat_f;
clear RGDP_hat_mnjf_france omega_tilde_mnjf_france p_hat_mnjf_france_rempat p_hat_mnjf_france

%% Rest of the countries

% calculating the cost of input bundle, b 
   
% Calculating P_hat for all countries
   
% P_hat_mnj: extending the matrix by J 
log_P_hat_mnj=log(P_hat_mnj);
log_P_hat_mnj0_repmat = kron(log_P_hat_mnj,ones(1,J));

% P_hat_weta = pi_M.*(P_hat_mnj_repmat)^(1-eta);
log_P_hat_weta = pi_M.*log_P_hat_mnj0_repmat;
  
% Summing by k and by i, so we just need to sum along the rows, this
% will give us a transposed vector [ 1xJ*N ] and the price of
% intermediate goods ==> P_M_hat
% Need to adjust for fact that some country-sector pairs never import,
% so will ahve a 0 P_M_hat raised to a negative number.
log_P_M_hat = squeeze(sum(log_P_hat_weta,1));
P_M_hat = exp(log_P_M_hat);   
 

% Multiplying by the intermediate share afer raising to 1-lambda
log_P_M_hat_wlambda = log_P_M_hat.*(1-pi_l)';  
   
% Repmat so that we have it by n, this means extending it by N
   
 log_P_M_hat_wlambda_repmat = repmat(log_P_M_hat_wlambda', [1 N]);

% Need to expand to mnj, but have to keep order so use a dummy (cant
% just use repmat here).
log_w_hat = log(wage_hat);
log_w_hat_to_vector=countrysecD*log_w_hat;
    
% Taking wage to power of 1-lambda and multiplying by labor share
   
log_w_hat_wlambda = pi_l.*log_w_hat_to_vector;

% Expanding the wage term by N so it is a [N*JxN]
% matrix
   
log_w_hat_wlambda_repmat = repmat(log_w_hat_wlambda,[1 N]);

% Create total cost: (N*J)*N matrix
log_b_hat_rep = (log_w_hat_wlambda_repmat + log_P_M_hat_wlambda_repmat);
b_hat_rep = exp(log_b_hat_rep);
b_hat = b_hat_rep(:,1);    

% p_hat_mnj: Calculating the denominator to calculate Q hat
p_hat_mnj = b_hat_rep.*repmat(a_hat,[1 N]);

% Q_hat_mnj: Taking change in exports divide by p_hat_mnj
Q_hat_mnj=X_hat_mnj./p_hat_mnj;
   
% VA_i_0_repmat: Replicating VA to expand it by N, since originally it
% is a [NJ 1] vector
VA_n_0_repmat = countrysecD*repmat(VA_n_0,[1 N]);
   
% omega_tilde_mnj_france: Weighting using pre shock VA
omega_tilde_mnj_france = VA_mnj_0./VA_n_0_repmat;

% Real GDP hat at the NJ x N level
RGDP_hat_mnj = omega_tilde_mnj_france.*Q_hat_mnj;
   
% Summing over importers
   
RGDP_hat_mj = nansum(RGDP_hat_mnj,2);
   
% RGDP_hat_m: real GDP for country m calculated using GDP deflator
% instead of CPI
RGDP_hat_m =countrysecD'*RGDP_hat_mj; 

% Clearing to save memory 
clear log_P_hat_mnj log_P_hat_mnj0_repmat log_P_hat_weta log_P_M_hat P_M_hat log_P_M_hat_wlambda;
clear log_P_M_hat_wlambda_repmat log_w_hat log_w_hat_to_vector log_w_hat_wlambda log_w_hat_wlambda_repmat;
clear log_b_hat_rep b_hat_rep b_hat p_hat_mnj;   
   
%% Assigning France to the RGDP

% RGDP_hat_m: Replacing real GDP for france with the number calculated
% using French firms
RGDP_hat_m(france) = RGDP_hat_m_france;

% GDP_deflator_i: Now we have the real gdp calculated using the GDP
% deflator we can back out the GDP deflator
GDP_deflator_m = VA_hat_n./RGDP_hat_m;

% RGDP_hat_f: Real GDP hat for each french firm.

RGDP_hat_f=VA_hat_f./GDP_deflator_m(france);

end

